Shifts in bird and butterfly communities of Tamhini Ghat

Introduction

This file was last update on 2021-06-03

Sampling effort

No. of visits to each transect/site.

abn%>%
  group_by(Study, Taxa, Transect)%>%
  summarise(n = length(unique(Date)))%>%
  spread(Transect, n)
Study Taxa T1A T1B T2 T3 T4 T5 T6
1998 - 2000 Birds 19 18 7 10 4 1 1
1998 - 2000 Butterflies 14 7 4 9 4 1 1
2016 - 2017 Birds 23 23 22 19 16 22 17
2016 - 2017 Butterflies 23 20 20 22 18 21 15

No. of individuals sampled

abn%>%
  group_by(Study, Taxa, Transect)%>%
  summarise(n = sum(Abundance))%>%
  spread(Transect, n)
Study Taxa T1A T1B T2 T3 T4 T5 T6
1998 - 2000 Birds 255 165 103 213 236 7 28
1998 - 2000 Butterflies 93 59 12 180 157 10 4
2016 - 2017 Birds 480 393 358 289 165 203 133
2016 - 2017 Butterflies 432 188 207 652 187 224 124

Uneveness in sampling effort across sites in studies may bias estimates of change.

Rarefying across sites

library(vegan)

# input for vegan

veg_birds <- abn%>%
  filter(Taxa == "Birds")%>%
  group_by(Study, Transect, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = paste(Study, "_", Transect))%>%
  group_by(sample)%>%
  #mutate(N = sum(n, na.rm = T))%>%
  #filter(N > 0)%>%#removing samples where no species were detected
  dplyr::select(sample, Specific.Epithet, n)%>%
  ungroup()%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

veg_butterfly <- abn%>%
  filter(Taxa == "Butterflies")%>%
  group_by(Study, Transect, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = paste(Study, "_", Transect))%>%
  group_by(sample)%>%
  #mutate(N = sum(n, na.rm = T))%>%
  #filter(N > 0)%>%#removing samples where no species were detected
  dplyr::select(sample, Specific.Epithet, n)%>%
  ungroup()%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")
# sample info

samp <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  dplyr::select(Study, Transect, Habitat)%>%
  distinct()%>%
  mutate(sample = paste(Study, "_", Transect))
## Species accumulation curves

rarebutts <- rarecurve(veg_butterfly) ## butterflies

# Coerce data into "long" form.
## taken from https://stat.ethz.ch/pipermail/r-sig-ecology/2018-December/005867.html

names(rarebutts) <- samp$sample

protox <- mapply(
  FUN = function(x, y) {
    mydf <- as.data.frame(x)
    colnames(mydf) <- "value"
    mydf$sites <- y
    mydf$subsample <- attr(x, "Subsample")
    mydf}, 
  x = rarebutts, y = as.list(names(rarebutts)), SIMPLIFY = FALSE)

rarebutts.df <- do.call(rbind, protox)

rownames(rarebutts.df) <- NULL  # pretty

rarebutts.df$Taxa <- "Butterflies"

## repeating for birds

rarebirds <- rarecurve(veg_birds)

# Coerce data into "long" form.
## taken from https://stat.ethz.ch/pipermail/r-sig-ecology/2018-December/005867.html

names(rarebirds) <- samp$sample

protox <- mapply(
  FUN = function(x, y) {
    mydf <- as.data.frame(x)
    colnames(mydf) <- "value"
    mydf$sites <- y
    mydf$subsample <- attr(x, "Subsample")
    mydf}, 
  x = rarebirds, y = as.list(names(rarebirds)), SIMPLIFY = FALSE)

rarebirds.df <- do.call(rbind, protox)

rownames(rarebirds.df) <- NULL  # pretty

rarebirds.df$Taxa <- "Birds"

## joining rarefaction values

rarefied <- full_join(rarebirds.df, rarebutts.df)

## adding sample info

rarefied <- left_join(rarefied, samp, by = c("sites" = "sample"))
# Plot.
ggplot(rarefied, aes(x = subsample, y = value, col = Study))+
  geom_line(size = 1)+
  labs(x = "Sample Size", y = "Species")+
  facet_grid(Taxa ~ Transect)

ggsave(last_plot(), filename = "./Figures and Tables/figureA.png", height = 8, width = 16, dpi = 300)  

Change in diversity across years

Species richness

Across the two studies

abn%>%
  group_by(Taxa, Study)%>%
  summarise(richness = length(unique(Specific.Epithet)),
            N = sum(Abundance, na.rm = T))
Taxa Study richness N
Birds 1998 - 2000 70 1007
Birds 2016 - 2017 105 2021
Butterflies 1998 - 2000 45 515
Butterflies 2016 - 2017 66 2014

More species of both taxa were encountered in current study compared to the previous study. Difference in sampling effort is apparent accross two studies.

Across sites in the two studies

No. of species encountered across sites sampled in the two studies.

abn%>%
  left_join(transects, by = c("Transect"))%>%
  group_by(Taxa, Study, Transect)%>%
  summarise(richness = length(unique(Specific.Epithet)))%>%
  spread(Transect, richness)
Taxa Study T1A T1B T2 T3 T4 T5 T6
Birds 1998 - 2000 47 37 28 36 29 4 11
Birds 2016 - 2017 66 59 53 47 31 38 28
Butterflies 1998 - 2000 27 22 8 37 24 5 3
Butterflies 2016 - 2017 41 32 36 44 31 39 27

Fewer species were encountered in sites with lower effort in the previous study, i.e., T5 and T6.

Shannon diversity

# Function to caluculate Shannon diversity

shanon <- function(c){
  
  if(sum(c) < 20){
    
    H <- "Sample size too small"
    
  }else{
    
    require(SpadeR)
    require(tidyr)
    
    div <- Diversity(c, datatype = "abundance")
    
    H <- data.frame(div$Shannon_diversity)%>%
      rownames_to_column("Method")
    
  }
  
  return(H)
  
}
# Calculating diversity in each site over years

div_year <- abn%>%
  group_by(Study, Taxa, Transect,  Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  spread(Specific.Epithet, n, fill = c(0))%>%
  group_by(Study, Taxa, Transect)%>%
  nest()%>%
  mutate(H = map(data, ~shanon(.)))%>%
  dplyr::select(-data)%>%
  unnest(H)%>%
  filter(grepl("MLE", Method) | is.na(Method))

# Summary

div_year%>%
  group_by(Taxa, Study)%>%
  skimr::skim(Estimate)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Study, mean, sd)

Variable type: numeric

Taxa Study mean sd
Birds 1998 - 2000 17.31 6.47
Birds 2016 - 2017 24.88 5.74
Butterflies 1998 - 2000 15.87 3.42
Butterflies 2016 - 2017 20.40 3.24
# t test

div_year%>%
  group_by(Taxa)%>%
  nest()%>%
  mutate(test = map(data, ~t.test(Estimate ~ Study, data = .)),
         sumry = map(test, broom::tidy))%>%
  dplyr::select(sumry)%>%
  unnest()%>%
  tstats()
Taxa 1998-2001 2016-2017 parameter statistic p.value
Birds 17.31433 24.87600 10.16882 -2.213178 0.0508557
Butterflies 15.86775 20.40014 6.07105 -2.152995 0.0742637
# plot

div_year%>%
  ggplot(aes(Taxa, Estimate, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5)+
  labs(y = "Effective number of species")+
  scale_fill_brewer(palette = "Greys")

ggsave(last_plot(), filename = "./Figures and Tables/figure2.png", height = 6, width = 8, dpi = 300)  

The first order diversity of both taxa has increased significantly across the two surveys.

Change in composition of species accross years

# sample info

samp_birds <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  filter(Taxa == "Birds")%>%
  dplyr::select(Study, Taxa, Transect, Habitat)%>%
  distinct()%>%
  mutate(sample = paste(Study, "_", Transect))

samp_butterfly <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  filter(Taxa == "Butterflies")%>%
  dplyr::select(Study, Taxa, Transect, Habitat)%>%
  distinct()%>%
  mutate(sample = paste(Study, "_", Transect))
# Similarity matrix by studies

sim_birds <- abn%>%
  filter(Taxa == "Birds")%>%
  group_by(Study, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = Study)%>%
  ungroup()%>%
  dplyr::select(sample, Specific.Epithet, n)%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

sim_butterfly <- abn%>%
  filter(Taxa == "Butterflies")%>%
  group_by(Study, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = Study)%>%
  ungroup()%>%
  dplyr::select(sample, Specific.Epithet, n)%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

# bray - curtis

data_frame(
  Taxa = c("Birds", "Butterflies"),
  Dissimilarity = c(as.numeric(vegdist(sim_birds, method = "bray")),
                    as.numeric(vegdist(sim_butterfly, method = "bray")))
)
Taxa Dissimilarity
Birds 0.5495376
Butterflies 0.6710162
# ordination

NMDS_birds <- metaMDS(veg_birds, distance = "bray", k = 2, model = "local")

NMDS_butterfly <- metaMDS(veg_butterfly, distance = "bray", trymax = 1000, k = 2, model = "local")


# Extracting scores and adding habitat vars

ord_birds <- data.frame(scores(NMDS_birds))%>%
  rownames_to_column("sample")%>%
  left_join(samp_birds, "sample")

ord_butterfly <- data.frame(scores(NMDS_butterfly))%>%
  rownames_to_column("sample")%>%
  left_join(samp_butterfly, "sample")

ord <- bind_rows(ord_birds, ord_butterfly)
#Testing change in composition of birds

adonis(veg_birds ~ Study, data = samp_birds)
## 
## Call:
## adonis(formula = veg_birds ~ Study, data = samp_birds) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Study      1   0.77547 0.77547  3.9632 0.24827  0.001 ***
## Residuals 12   2.34798 0.19567         0.75173           
## Total     13   3.12345                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Testing change in composition of butterflies

adonis(veg_butterfly ~ Study, data = samp_butterfly)
## 
## Call:
## adonis(formula = veg_butterfly ~ Study, data = samp_butterfly) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Study      1    0.9191 0.91908   4.054 0.25253  0.001 ***
## Residuals 12    2.7205 0.22671         0.74747           
## Total     13    3.6396                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# plotting

ord%>%
  ggplot(aes(NMDS1, NMDS2, col = Study, fill = Study))+
  stat_ellipse(level = 0.95, geom = "polygon", linetype = "dashed", size = 1, alpha = 0.3)+
  geom_point(aes(shape = Habitat), size = 3)+
  guides(shape = guide_legend(nrow = 2, byrow = T))+
  scale_color_brewer(palette = "Set1")+
  scale_fill_brewer(palette = "Set1")+
  facet_wrap(~Taxa)+
  theme(legend.box = "vertical")

ggsave(last_plot(), filename = "./Figures and Tables/figure3.png", height = 6, width = 8, dpi = 300)

Both bird and butterfly communities compositions are significantly different compared to the previous study.

Change in niche width of community

Trophic Guilds

# compiling trophic guild data

host <- read.csv("./Data/Butterflies_host plant_280820.csv") #butterfly host plant data

diet <- read.csv("./Data/Birds_diet.csv") #bird diet data

trophic <- host%>%
  full_join(diet)%>%
  dplyr::select(Specific.Epithet, Guild)%>%
  distinct()

# Calculating diversity of guilds

guilds <- abn%>%
  left_join(trophic, by = c("Specific.Epithet"))%>%
  filter(!Guild == "")%>%
  group_by(Study, Taxa, Transect, Guild,  Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  spread(Specific.Epithet, n, fill = c(0))%>%
  group_by(Study, Taxa, Transect, Guild)%>%
  nest()%>%
  mutate(H = map(data, ~shanon(.)))%>%
  dplyr::select(-data)%>%
  unnest(H)%>%
  filter(grepl("MLE", Method) | is.na(Method))

# summary

guilds%>%
  group_by(Taxa, Guild, Study)%>%
  skimr::skim(Estimate)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Guild, Study, mean, sd)

Variable type: numeric

Taxa Guild Study mean sd
Birds Carnivore 1998 - 2000 3.43 0.31
Birds Carnivore 2016 - 2017 6.26 0.15
Birds Frugivore 1998 - 2000 4.11 0.92
Birds Frugivore 2016 - 2017 4.21 0.60
Birds Grainivore 1998 - 2000 3.86 0.75
Birds Grainivore 2016 - 2017 3.84 0.49
Birds Insectivore 1998 - 2000 8.77 1.95
Birds Insectivore 2016 - 2017 13.58 4.82
Birds Omnivore 1998 - 2000 4.82 0.92
Birds Omnivore 2016 - 2017 5.88 1.15
Butterflies Generalist 1998 - 2000 3.16 1.48
Butterflies Generalist 2016 - 2017 5.95 1.87
Butterflies Grass Specialist 1998 - 2000 1.80 NA
Butterflies Grass Specialist 2016 - 2017 2.71 0.71
Butterflies Herb Specialist 1998 - 2000 3.58 NA
Butterflies Herb Specialist 2016 - 2017 4.09 0.86
Butterflies Liana Specialist 1998 - 2000 NaN NA
Butterflies Shrub Specialist 1998 - 2000 4.36 0.77
Butterflies Shrub Specialist 2016 - 2017 4.18 2.10
Butterflies Tree Specialist 1998 - 2000 5.92 3.07
Butterflies Tree Specialist 2016 - 2017 6.88 1.63
# t - test

guilds%>%
  filter(Guild != "Liana Specialist")%>% #removing due to lack of obs
  complete(nesting(Study, Taxa), Transect, Guild, fill = list(H = 0))%>%
  group_by(Taxa, Guild)%>%
  nest()%>%
  mutate(test = map(data, ~t.test(Estimate ~ Study, data = .)),
         sumry = map(test, broom::tidy))%>%
  dplyr::select(sumry)%>%
  unnest()%>%
  tstats()
Taxa Guild 1998-2001 2016-2017 parameter statistic p.value
Birds Carnivore 3.427333 6.257500 28.84235 -43.1080663 0.0000000
Birds Frugivore 4.114600 4.211286 55.19450 -0.5971850 0.5528259
Birds Grainivore 3.856667 3.840000 32.48672 0.0998621 0.9210676
Birds Insectivore 8.766000 13.578714 68.06664 -6.6782693 0.0000000
Birds Omnivore 4.821000 5.877000 81.46157 -5.0537422 0.0000026
Butterflies Generalist 3.162500 5.946333 21.02601 -6.3513785 0.0000027
Butterflies Grass Specialist 1.798000 2.710400 34.00000 -8.3567205 0.0000000
Butterflies Herb Specialist 3.581000 4.092143 48.00000 -4.4413461 0.0000524
Butterflies Shrub Specialist 4.357500 4.179400 44.86704 0.4999673 0.6195415
Butterflies Tree Specialist 5.922000 6.879714 26.24765 -1.5944504 0.1228073
# plot
guilds%>%
  ggplot(aes(reorder(Guild, Estimate), Estimate, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5, position = position_dodge(preserve = "single"))+
  labs(y = "Effective number of species", x = "Trophic Guild")+
  scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  scale_fill_brewer(palette = "Greys")+
  facet_wrap(~Taxa, scales = "free_x", ncol = 1)

ggsave(last_plot(), filename = "./Figures and Tables/figure4.png", height = 9, width = 8, dpi = 300)  

Carnivorous, Grainivorous birds have reduced in the comunity, compared to insectivorous birds which have increased. Similarly, Shrub specialist and tree specialist butterflies have reduced where grass specialists and herb specialist have increased. Suprsingly, generalist species of both taxa appear to be unchanged.

Change in niche specialisation

-Can we incorporate abunndance into community specialisation index?

Trophic niche

# trophic specialisation based on Juliard et al. 2006

butterflies_TSI <- host%>%
  mutate(H = length(unique(Hostplant.Family)))%>%
  group_by(Specific.Epithet)%>%
  mutate(h = length(unique(Hostplant.Family)))%>%
  ungroup()%>%
  mutate(TSI = sqrt((H/h)-1))

birds_TSI <- diet%>%
  mutate(H = length(unique(Prey)))%>%
  group_by(Specific.Epithet)%>%
  mutate(h = length(unique(Prey)))%>%
  ungroup()%>%
  mutate(TSI = sqrt((H/h)-1))

CSI.T <- abn%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet)%>%
  group_by(Study, Taxa, Transect)%>%
  distinct()%>%
  left_join(butterflies_TSI, by = c("Specific.Epithet"))%>%
  left_join(birds_TSI, by = c("Specific.Epithet"))%>%
  mutate(TSI = ifelse(is.na(TSI.x), TSI.y, TSI.x))%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet, TSI)%>%
  summarise(CSI.T = mean(TSI, na.rm = T))
# summary
CSI.T%>%
  group_by(Taxa, Study)%>%
  skimr::skim(CSI.T)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Study, mean, sd)

Variable type: numeric

Taxa Study mean sd
Birds 1998 - 2000 2.59 0.16
Birds 2016 - 2017 2.53 0.09
Butterflies 1998 - 2000 5.34 0.39
Butterflies 2016 - 2017 5.07 0.24
# t - test

CSI.T%>%
  group_by(Taxa)%>%
  nest()%>%
  mutate(test = map(data, ~t.test(CSI.T ~ Study, data = .)),
         sumry = map(test, broom::tidy))%>%
  dplyr::select(sumry)%>%
  unnest()%>%
  tstats()
Taxa 1998-2001 2016-2017 parameter statistic p.value
Birds 2.588009 2.531490 9.462567 0.8103492 0.4376441
Butterflies 5.338563 5.065191 10.056310 1.5727722 0.1466765
CSI.T%>%
  ggplot(aes(Taxa, CSI.T, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5)+
  labs(y = "Community specialisation index \n (Trophic)")+
  scale_fill_brewer(palette = "Greys")+
  scale_y_sqrt()

Niche width of the community did not change over the two studies.

Habitat Niche

HSI <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  group_by(Study, Taxa, Habitat, Specific.Epithet)%>%
  # calculating abundance of each species in each habitat type
  summarise(n = sum(Abundance, na.rm = T))%>%
  group_by(Study, Taxa)%>%
  mutate(N = sum(n), # total number of indviduals encountered 
         rel.prop = n/N)%>% # relative proportion of species
  group_by(Study, Specific.Epithet)%>%
  summarise(HSI = sd(rel.prop)/mean(rel.prop))%>%
  replace_na(replace = list(HSI = 0))

CSI.H <- abn%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet)%>%
  group_by(Study, Taxa, Transect)%>%
  distinct()%>%
  left_join(HSI, by = c("Study", "Specific.Epithet"))%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet, HSI)%>%
  summarise(CSI.H = mean(HSI, na.rm = T))
# summary
CSI.H%>%
  group_by(Taxa, Study)%>%
  skimr::skim(CSI.H)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Study, mean, sd)

Variable type: numeric

Taxa Study mean sd
Birds 1998 - 2000 0.54 0.10
Birds 2016 - 2017 0.49 0.02
Butterflies 1998 - 2000 0.61 0.09
Butterflies 2016 - 2017 0.57 0.03
# t - test

CSI.H%>%
  group_by(Taxa)%>%
  nest()%>%
  mutate(test = map(data, ~t.test(CSI.H ~ Study, data = .)),
         sumry = map(test, broom::tidy))%>%
  dplyr::select(sumry)%>%
  unnest()%>%
  tstats()
Taxa 1998-2001 2016-2017 parameter statistic p.value
Birds 0.5363977 0.4942803 6.719090 1.068760 0.3220632
Butterflies 0.6110527 0.5692785 7.782004 1.210166 0.2616879
CSI.H%>%
  ggplot(aes(Taxa, CSI.H, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5)+
  labs(y = "Community specialisation index \n (Habitat)")+
  scale_fill_brewer(palette = "Greys")

Habitat niche of the community did not change over the study periods.

Seasonal changes in relative abundance

season <- abn%>%
  mutate(Month = month(dmy(Date), label = T))%>%
  group_by(Study, Month, Taxa)%>%
  summarise(n = sum(Abundance, na.rm=T))%>%
  group_by(Study, Taxa)%>%
  mutate(N = sum(n),
         rel.prop = n/N)

season%>%
  ggplot(aes(Month, rel.prop, col = Taxa, shape = Taxa))+
  geom_point(size = 3)+
  geom_path(aes(group = Taxa), size = 1, linetype = "dashed")+
  scale_color_brewer(palette = "Set1")+
  facet_wrap(~Study, ncol = 1)

Site Map

library(raster)
library(sf)

# Raw point data

points <- read.csv("./Data/sites.csv")%>%
  drop_na()%>%
  dplyr::select(-Location)

# Converting to SF points

points_sf <- st_as_sf(points, coords = c('Lon', 'Lat'))

st_crs(points_sf) <- 4326

# Making transect lines

lines <- points_sf%>%
  mutate(Transect = as.factor(Transect))%>%
  group_by(Transect)%>%
  summarise()%>%
  st_cast("LINESTRING")
# getting study extent

study_ext <- st_bbox(points_sf)

ext <- st_as_sfc(study_ext)

# getting osm map

library(osmdata)

# Tamhini WLS

tamhini_raw <- opq(bbox = study_ext, timeout = 100)%>%
  add_osm_feature(key = 'boundary', value = 'protected_area')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
tamhini <- tamhini_raw$osm_multipolygons

st_crs(tamhini) <- 4326

# Roads

tamhini_area <- st_bbox(tamhini)

roads_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
  add_osm_feature(key = 'highway')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
roads <- roads_raw$osm_lines

st_crs(roads) <- 4326

# Landuse

landuse_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
  add_osm_feature(key = 'natural')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
forest <- landuse_raw$osm_polygons%>%
  filter(natural == "wood")

water <- landuse_raw$osm_multipolygons

water_small <- landuse_raw$osm_polygons%>%
  filter(natural == "water")

st_crs(forest) <- 4326

# waterways

waterway_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
  add_osm_feature(key = 'waterway')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
waterways <- waterway_raw$osm_lines

st_crs(waterways) <- 4326

# elevation 

library(elevatr)

tam_elev <- get_elev_raster(locations = tamhini, z = 8) # fix contours
## Error in if ((is.null(locs) | is.na(locs)) & is.na(prj)) {: argument is of length zero
tam_contours <- rasterToContour(tam_elev, levels = seq(0, 1000, 100))
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'sampleRegular': object 'tam_elev' not found
# making final map

library(gridExtra)

fig1 <- grid.arrange(fig1a, fig1b, fig1c, legend, layout_matrix = rbind(c(1, 1, 1),
                                                          c(1, 1, 1),
                                                          c(2, 3, 4)))

ggsave(fig1, filename = "./Figures and Tables/figure1.png", height = 8, width = 8)